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Concentration gradients of signaling molecules are essential for patterning during development 
and they have been observed in both unicellular and multicellular systems. In subcellular systems, 
clustering of the signaling molecule has been observed. We develop a theoretical model of cluster- 
mediated concentration gradient formation based on the Becker-Doring equations of aggregation- 
fragmentation processes. We show that such a mechanism produces robust concentration gradients 
on realistic time and spatial scales so long as the process of clustering does not signihcantly stabilize 
the signaling molecule. Finally, we demonstrate that such a model is applicable to the pomlp 
subcellular gradient in fission yeast. 

PACS numbers: 87.10.-e, 87.16.A-,87.17.Aa,87.18.Tt 


I. INTRODUCTION 

Individual cells make precise cell-fate decisions based 
on information from signaling networks during develop¬ 
ment. A central paradigm of information transfer during 
development is the morphogen gradient, a spatially vary¬ 
ing concentration profile [T] . The dynamics of morphogen 
gradient formation [2] and interpretation [3] has been the 
subject of intense study, both experimentally mi] and 
theoretically [SHU]- 

Single cells can interpret spatial gradients across the 
scale of the cell itself, either in response to external chem¬ 
ical gradients - such as chemotactic gradients m - or 
subcellular gradients. Subcellular concentration gradi¬ 
ents are found in the single cell stage of C. elegans devel¬ 
opment [HI, and in single cell organisms including fission 
yeast [m [H] and bacteria (HMI]- In parallel with ex¬ 
perimental work, theoretical models have demonstrated 
how subcellular gradients can be formed on the relevant 
time and spatial scales [22H1I1- 

Clustering, either of the signaling molecules or recep¬ 
tors, is observed in a number of subcellular signaling 
systems PHMT] . Receptor clustering helps ensure reli¬ 
able readout of input signal, for example through recep¬ 
tor clustering by positive feedback that enables binary 
on/off decisions [35]. Clustering of signaling molecules is 
also observed [55], potentially altering the dynamics of 
concentration gradient formation. 

Here, we develop a mechanistic model of cluster- 
mediated concentration gradient formation based on the 
Becker-Doring equations of aggregation-fragmentation 
processes [SS] [31] . We show that such a model can pro¬ 
duce a robust concentration gradient under (certain) bi¬ 
ological relevant parameter conditions. In particular, we 
find that the process of clustering must not significantly 
increase the signal molecule effective lifetime. Finally, 
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we apply the clustering model to the cortical subcellu¬ 
lar gradient pomlp in fission yeast [HI [M [551 [33]. The 
modelling describes how a single component (and its in¬ 
teractions with itself) can create a robust subcellular 
concentration gradient by adapting its effective dynami¬ 
cal properties at different spatial positions. Theoretical 
approaches on the dynamics of concentration gradient 
formation in subcellular systems may need to be signifi¬ 
cantly different from embryonic systems. 


II. AGGREGATION-FRAGMENTATION 
MODEL OF GRADIENT FORMATION 

We develop a Becker-Doring-like model of concentra¬ 
tion gradient formation via clustering. This model takes 
into account aggregation, fragmentation and diffusion of 
clusters with only a single molecular species. The param¬ 
eters used are derived from experiments in fission yeast 
[29] ■ with typical length and time scales on the order of 
a few microns and seconds respectively. 


A. Model motivation 

We consider the Becker-Doring equations with con¬ 
served number and diffusion |36II38j . The motivation 
for such a formalism comes from studies of pomlp in 
fission yeast. In time lapse movies, clusters of pomlp 
are not observed to coalesce (as in a Smoluchowski pro¬ 
cess [34] , which describes E-Cadherin clustering [35] ) but 
observed to grow and decay on second time scales, ap¬ 
parently independently of other clusters [55] ■ Therefore, 
we assume that only monomeric molecules are taken up 
into a cluster and during fragmentation single molecules 
are released from clusters, Fig. Cluster disassocia- 
tion events from the membrane are not observed in cells 
[55] and so we assume that disassociation involves only a 
single molecule at a time (ie. when a component disas¬ 
sociates from both the cluster and membrane it does so 
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without affecting the other components of the cluster), 
Finally, large clusters are typically not observed 
to join at the insertion region so we consider inser¬ 
tion only in the monomeric form. While difficult to solve 
analytically the existence and uniqueness of solutions of 
qualitatively similar models can be shown in relevant pa¬ 
rameter regimes 1371. We use a one-dimensional model of 
concentration gradient formation as our aim is to high¬ 
light the important general behavior of such a model. 


B. Model formulation 

We consider a mean-field discrete model of cluster for¬ 
mation. The reaction schemes are shown in Fig. and 
result in the following equations for clusters where 
represents the concentration of a cluster containing s 
molecules: 

drii d'^rii 

Srnax 

- ^ js - j2 + JSix) ( 1 ) 

dus d^Us . 

- {js+i - ^J,s+l{s + l)ns+i) for s > 1 (2) 

js =as-inins-i - PsTIs ■ ( 3 ) 

Ds denotes the diffusion coefficient, fXs the membrane dis- 
association rate and as, Ps denote the aggregation and 
fragmentation rates respectively for a cluster containing s 
molecules. J is the monomeric insertion rate and bound¬ 
ary conditions are ^|^|a;=o.z, = 0 for s > 2. 


C. Parameter scaling 
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FIG. 1. (Colour online) A Reaction scheme of processes in 
the clustering mechanism, Eqs. m Crossed circles denote 
particles exiting the system {e.g. through membrane disas¬ 
sociation). B Scaling of dynamic parameters as function of 
cluster size. For aggregation the dashed line corresponds to 
the scaling of a without the term . For disassociation, 

the black line corresponds to equal rate of disassociation of 
particles in all cluster sizes. Gray line corresponds to = 1. 


Eqs. [T][^ have a large parameter space, increasing with 
Smax, tbe maximum cluster size. However, using bio¬ 
physical arguments we can reduce the parameter space 
to seven non-dimensional parameters as discussed below. 

Diffusion: We take Ds = Dis~^^ where = 1, 
Fig. [Tp, consistent with experiments on clustered pro¬ 
tein membrane diffusion [331110] (note, this contrasts with 
theoretical predictions of logarithmic diffusion scaling on 
biological membranes |41j ). We have confirmed our main 
conclusions hold for rjj = 2/3 (not shown). 

Disassociation: Experimentally, large clusters are not 
observed to disassociate from the membrane[33]- There¬ 
fore, we assume only single molecules disassociate in each 
reaction. We take /is = /iis“’'^‘. = 0 if clustering has 

no membrane-stabilizing effect {i.e. the rate of monomer 
disassociation is equivalent from all clusters). We also 
consider the case = 1, i.e. clustering stabilizes the 
signaling molecule within the cluster. Fig. We shall 
see that the latter scenario results in a non-robust gradi¬ 
ent and is a key result of this analysis. 


Aggregation: Aggregation depends on the cluster size 
- bigger clusters are more likely to collide with and ag¬ 
gregate a monomer. How exactly aggregation depends 
on the cluster topology is unclear and here we consider 
Us = with Ta = 1/2. Experimentally, there is an 

upper limit on the cluster size |29j , so an additional term 
within the aggregation parameter is included to limit the 
maximum size of the clusters: Og = where 

So = 100 [33], Fig- Such limitation on the maximum 
cluster size could also have been implemented by includ¬ 
ing a component in the fragmentation rate (below). 

Fragmentation: Scaling of the fragmentation rate is 
likely to occur via a similar mechanism to aggregation. 
Hence, we consider Ps = Pis''^ with rp = ra = 1/2. 
We have confirmed that using „ = 1/3 or rp^ = 2/3 
does not significantly alter our key conclusions so long 
as rp = ra. We choose the values of ai and Pi such 
that cluster dynamics are on the order of a few seconds, 
consistent with pomlp cluster dynamics [33] ■ Note, we 
use Pi as the lowest order for fragmentation (even though 
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only terms Ps, s > 2 appear in Eqs. T|3) for clarity in 
representing the scaling. 

Non-dimensionalizing by substituing t = TjfXi, Ps = 
^Dini ^ “ UyJDi/pi, along with the scaling ar¬ 

guments, reduces Eqs.flll^to 


dtpi _ d‘^(j)i 
dr du^ ^ 


-91 + 


s=2 

ol-r 






’^*02 +/3</>2 — K</>1 +(5(u) (4) 


9(/)s 


- S -I- (s -I- 1) 


l-J-u 


s-l-1 


dr du^ 

- (s - 

-/3 - (s + lY'^(l)s+l) , 


(5) 


where n = Ja\/Dp^ and /3 = Pijpi. The large param¬ 
eter space in Eqs. m has been reduced to seven dimen¬ 
sionless parameters (Table |T|, independent of Smax- Four 
of these seven parameters, ra,p,D,^l are constrained by 
physical arguments as described above. The phenomeno¬ 
logical cluster size factor sg is limited by experimental 
observation of cluster sizes [29]. Only k and /? are free 
parameters and hence, despite the apparent complexity 
of Eqs. E13 the dynamic behavior of the system is effec¬ 
tively described by just two parameters. 


D. Concentration gradient profile 


We solve the above system of equations (up to s^ax 
cluster size, typically 200) using Matlab pde45 and con¬ 
firm the steady-state distribution using the Matlab ode 
solver bvp4c. The total concentration profile, defined as 


TVt = ^ sus , (6) 


is shown in steady-state in Fig. UK, for a range of k. 
Unsurprisingly, large clusters are localized to the source 
region, Fig. [^. The effective diffusion. 
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FIG. 2. (Colour online) A Concentration profile, Eq.j^with 
K = 425 and /3 = 20. Dotted and dashed black lines cor¬ 
respond to K « 142 (J —>■ J/3) and k — 1275 (J —>■ 3J) 
respectively. Gray lines show solutions of a model without 
clustering with same insertion and disassociation rates and 
D = 0.1pms~^ (see SDD model in Table |n|. Inset shows pro¬ 
files on logarithmic scale. B Cluster concentration as function 
of position, with darker color representing higher concentra¬ 
tion. Red (step-like) and blue (smooth) lines correspond to 
the cluster size with highest concentration and mean cluster 
concentration respectively. Inset: effective diffusion, Eq. as 
function of position for j3 = 20 (black), /3 = 200 (dashed) and 
P = 2 (gray). Parameters as Fig.unless otherwise stated. 


TABLE I. Dimensionless parameters in clustering model 


Parameter 

Value 

Note 

^ a. 1 So 

Y 100 

So = 100 constrain 
maximum cluster size 

rp 

ri^ 

1 

2 

0,1 

r/3 = Va 

Tfj, > 0 results in stabilization 
of monomers in larger clusters 

rn 

1 

See [291 HO] 

„ — Jatl 

P 

~ 425 

~ 20 

Defines range of clustering effects 

Fragmentation occurs on shorter 



timescales than disassociation 


is a function of position, increasing away from the source. 
Fig. inset. The timescales of clustering [42] and pro¬ 
tein dynamics are consistent with observations of pomlp 
in vivo [29] . 

What biophysical processes do k and /3 represent? f3 
effectively defines the relative lifetime of molecules in 
clusters before disassociation, which in turn alters the 
effective diffusion. Fig. inset, k compares the effects 
of insertion and aggregation (increasing either amplifies 
clustering) with diffusion and disassociation (increasing 
either reduces clustering). Small k corresponds to weak 
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clustering and the diffusive dynamics of the monomer 
dominate. Large k results in clustering dominating the 
dynamics with a resulting steep concentration gradient. 
Using biologically plausible parameter values [22]) both 
K and /3 occur at values that allow large clusters to form 
(to give benefits of modulating diffusion) without per¬ 
mitting very large clusters to dominate the dynamics, 
particularly away from the source region. 


III. ROBUSTNESS OF CLUSTER-MEDIATED 
CONCENTRATION GRADIENT FORMATION 

We have demonstrated that a clustering model can pro¬ 
duce a concentration gradient on similar spatial and tem¬ 
poral scales observed in subcellular systems. Ultimately, 
the concentration profile must be able to impart precise 
(j.e. robust) information to the cell. Below, we explore 
whether such concentration gradients can be robust to 
relevant biological fluctuations. To gain a qualitative un¬ 
derstanding, we first discuss a phenomenological model 
incorporating concentration-dependent diffusion before 
discussing the effects of noise on the clustering model. 


A. Concentration-dependent diffusion 

We consider a one-dimensional reaction-diffusion equa¬ 
tion with concentration-dependent diffusion (CDD), 
D{p), a function of the local protein concentration, p\ 



with boundary conditions D{p(t)))^\^^Q -p J = 0 where 
J is the protein insertion rate and p{x —>■ oo) = 0. We 
consider the case D[p) = Do{p/p)'^. The solution to 
Eq. [^then has steady-state solution p{x) = B{x + xi) 
where m = ^jr^ B = A’"p[m(m — 1)]™/^ and Xi = 

1 T71-2 

A j ^ 5(”T-)) where g{m) = (m — 1)) ^ 

Po = J/'JDqp and A = Table [u} 

In Fig. iA we compare steady-state profiles for the 
NLD and CDD models (right inset shows the behavior of 
the diffusion coefficient in the CDD model as a function of 
position) with the scenario of linear diffusion {D{p) = Dq 


TABLE II. Steady state solutions to SDD, NLD and CDD 
models. For NLD model solution for non-linear degradation 
term —ap^ is shown (9). CDD model solution for r = 1/2. 


Model 

p{x, t — >■ oo) 


Parameters 


SDD 

J 

V-^0^^0 


A = \/ — ~ 1.6/im 

V MO ^ 


NLD 

A 

(a:+a;o)^ 

A = 

6Do/a, xo = 

^1/3 

CDD 

B 

{x+xi)‘^ 

B = 

144M^ *1 = A 

jl/3 



FIG. 3. A Steady-state profiles for SDD (gray), NLD 
(dashed) and CDD (black) models (Table |^. SDD model: 
J = 16.7s“'^, Do = and p = 0.03s“ . p = 40pm~^ 

and non-linear degradation rate a were chosen such that the 
total particle number was equal in all models with same J, 
Do and po and such that the average diffusion in the CDD 
model for 0 < a;/A < 3 equals Do- Left inset: same as main 
panel but on y-log scale. Right inset: D{p) in CDD model, 
dashed line corresponds to D = O.lpm^s~^. B Error in defin¬ 
ing spatial position, Eq. due to variations in J, given by 
SJ = O.IJ. Black bars denote parameters as A. Clear bars 
correspond to equal characteristic length scale, X = xo = 

{p — 6pm~^). Smaller values denote more robust profiles. C 
Robustness of the NLD and CDD models due to intrinsic fluc¬ 
tuations at a:/A = 2 relative to SDD model. Black and clear 
bar notation same as B. Though the specific fractions vary as 
a function of position, the trend of Sxf^^ < Sx'pJ^^ < 
typically holds in the region 0 < a; < 3A. 


in Eq. and degradation (SDD model). Table [TT] We 
define robustness as the positional error, 5x in defining a 
boundary at a threshold concentration [^ 

5x{x, t) = Sp{x, t)/\p{x, f)'| (9) 

due to concentration fluctuations 5p{x,t) {p' is the spa¬ 
tial derivative) from, for example, variations in J. Sx is 
dependent on when and where measurement occurs. We 
focus on the spatial position in steady-state for distances 
around 2 — 3pm from the source, consistent with bound¬ 
aries in fission yeast. Interestingly, in steady-state, dx is 
independent of position for all three models if the concen¬ 
tration fluctuations are due to variation in the injection 
rate J, Fig. 3B, [Sj. Previously, it has been shown that 
a model with non-linear degradation (NLD) can produce 
concentration profiles that are robust to variations in the 
insertion rate |2]. In Fig. we demonstrate that the 
CDD model, just as with the NLD model, is more robust 
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to variations in the insertion rate compared to the SDD 
model when 0 < r < 1. 

Another source of error are stochastic biochemical (in¬ 
trinsic) fluctuations, relevant in both embryonic and sub- 
cellular systems [251113H1S] • Such fluctuations are typi¬ 
cally well-described by Poisson statistics [46l |47] : 


Spint — Oy jjrp ) (10) 

where a is a constant that is assumed to be model inde¬ 
pendent and T is the averaging period. NLD models are 
generally less robust to intrinsic fluctuations m- How¬ 
ever, in the CDD model the diffusion coefficient increases 
with distance from the source, which in turn increases 
the effects of time averaging [46] and hence reduces the 
detrimental effects of intrinsic fluctuations. In Fig.[^ we 
show that the SDD model is most robust to such varia¬ 
tions at x/\ = 2, but concentration-dependent diffusion 
results in more robust gradients than those formed by 
non-linear degradation processes. 


B. Robustness of clustering model 

Having developed a qualitative understanding of the 
robustness of concentration-dependent diffusion to rele¬ 
vant fluctuations we now discuss the robustness of the 
full clustering model, Eqs[g^ 

To test the robustness of the clustering model to varia¬ 
tions in protein insertion we first created 250 profiles us¬ 
ing parameters in Fig.[T]but with the insertion rate Gaus¬ 
sian distributed (mean J, standard deviation 5J = 0.2J), 
Fig. Fig. shows that the concentration at x = 0 
(normalized by the value when J = J) increases rapidly 
(faster than equivalent SDD model) as J increases; as 
expected, the fluctuations at the source are increased in 
the clustering model. Following |9], we define the char¬ 
acteristic length scale of the profile as 


the source increase especially quickly with increasing 5J. 
Here, our more precise analysis shows that clustering re¬ 
sults in a robust gradient when comparing the competing 
effects of p{x = 0) and Ae// fluctuations with changing 
J (note that the particular effect is position dependent). 
Furthermore, we see that the robustness is due to the 
effective spatial adaptation of the diffusion coefficient, 
which is highly dependent on ,5, Fig. [^, as well as on 
the scaling of the disassociation rate with cluster size. 
The simplifications used in the two-state model meant 
that the critical relationship between fragmentation and 
disassociation rates was not appreciated, highlighting a 
further advantage of our more detailed approach. 

Finally, we compare the effects of intrinsic fluctuations 
on the robustness of the clustering model. The processes 
of clustering are non-linear and therefore the Poisson ap¬ 
proximation is less valid. However, direct experimen¬ 
tal measurement of intrinsic fluctuations in a subcellular 
gradient [29] suggest that away from the source region 
{x > 2^m), where diffusion is the predominant dynamic 
mechanism, that intrinsic fluctuations are approximated 
closely by Poisson statistics. Therefore, Eq. 10 is a rea¬ 


sonable approximation to the intrinsic noise here, though 
likely represents a lower bound on the true intrinsic noise 
due to neglected non-linear effects from clustering pro¬ 
cesses. In Fig. EP we show the positional accuracy of the 
cluster model given Poisson distributed intrinsic fluctu¬ 
ations, compared to the SDD model with similar profile 
shape and protein disassociation rate. The two models 
have qualitatively similar sensitivity to intrinsic fluctua¬ 
tions due to their similar profile shape, Eq. [T^ The clus¬ 
tering model is less sensitive to intrinsic fluctuations than 
the CDD and NLD models discussed above as the latter 
two models have algebraic, rather than exponential-like, 
profiles [47] . 


IV. POMIP SUBCELLULAR GRADIENT 


= pix)/\p'{x)\. (11) 

Ae// is a function of position and henceforth we consider 
the mean Ag// in the range 2/xm < x < 3/rm, 

In Fig. E^, we see that (A) 2 - 3 ^m scales inversely with 
increasing J (unlike the SDD model). Therefore, the 
increases variation at x = 0 is compensated for by adap¬ 
tation in the characteristic profile length. This results 
in the cluster model providing precise spatial informa¬ 
tion, Fig. 1^, compared to SDD model. In Fig. [^ we 
also show that the positional accuracy of the clustering 
model is significantly reduced when clustering stabilizes 
the protein. Therefore, the clustering model is only ro¬ 
bust if individual monomers only spend a relatively short 
period in each cluster - if they are too stable then the ro¬ 
bustness is lost. 

The trend of Ae// being inversely proportional to J 
can be derived using a two-state (monomer / clustered) 
model [25]. However, in such a model fluctuations at 


We apply our clustering model to the specihc case of 
pomlp in fission yeast and its repression of the down¬ 
stream target cdr2p [T5l[T6l|48]. We incorporate spatially 
distributed pomlp insertion in the polar region mm- 
The formation of the pomlp subcellular gradient has 
been modelled previously [la Hi EO]. However, these 
approaches either did not consider clustering [451 EO] or 
presented only a qualitative model of clustering with only 
two states [25] . 

Cdr2p itself is known to cluster; indeed, it forms sig¬ 
nificantly larger clusters with around 80-100 molecules 
in the clusters localized to the cell center m- Pomlp 
has a dual effect on cdr2p. First, it is involved in cdr2p 
dephosphorylation, resulting in cdr2p membrane disasso¬ 
ciation. Second, it inhibits cdr2p cluster formation. This 
double mechanism of repression helps produce the sharp 
response of cdr2p to pomlp inhibition [48], Fig. [^,B. 
We use a phenomenological model of cdr2p cluster for¬ 
mation since our focus is on pomlp and how accurately it 
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FIG. 4. A 250 profiles (gray) generated from parameters 
shown in Fig. , but with 20% fluctuations in insertion rate. 
Black line is mean concentration. Inset is same with y-axis 
on logarithmic scale. B Scaling of concentration at a; = 0 
(triangles for clustering model, dashed line for SDD model) 
and decay length {X) 2 - 3 tim (squares for clustering model, line 
for SDD model), see text for definition, as function of inser¬ 
tion rate normalized by their value for mean insertion rate, 
denoted by J, po and A. C Positional error of the cluster¬ 
ing model relative to the SDD model with equivalent kinetic 
parameters. Black line corresponds to clustering model with 
equal disassociation rate in all clusters, = 0. Gray line cor¬ 
responds to clustering model with = 1. D The positional 
accuracy of the clustering model in presence of intrinsic fluc¬ 
tuations, Eq.|10[ normalized by the equivalent accuracy of the 
SDD model at x = 0. Equivalent result for SDD model shown 
as gray line. 

can define the cdr2p boundary. Therefore, we use a two 
state model for cdr2p (monomeric or clustered, similar to 
[29| ) with additional interactions between pomlp (where 
[P] denotes the total pomlp concentration at a particular 
position, regardless of the particular cluster distribution, 
Eq. 1^ and cdr2p (where [Ci^ 2 ] denotes cdr2p concen¬ 
tration in monomeric and clustered forms respectively), 
Fig. [^. Cdr2p is assumed to be inserted uniformly (rate 
Jc) across the membrane m- 

+ Iic. 2 [C 2 ] - [P](7l[Cl] - 72[C2]) + Jc (12) 

d[C2] ^ d^[C2] rri , , qpiurG 

—Dc,2—^ - ^J-c,2[C2\ + ac([PJ)[GiJ 

- I3c[C2]-12[P][C2] (13) 

The direct inhibition of cdr2p on the cortex by pomlp 
is approximated by the term — 7 i, 2 [C'i, 2 ][P] and the ag¬ 
gregation factor Qfc is now dependent on the concentra¬ 
tion of pomlp. Fig. [^. We assume a Hill-like behav- 


FIG. 5. (Colour online) A Experimentally measured (nor¬ 
malized) intensity profiles for pomlp (red) and cdr 2 (green) 
on the cortex of fission yeast cells (n = 14), see ^T] for de¬ 
tails. Distance defined as the distance around the cortex of 
the cell starting from the tip, as shown by the white line in 
inset. Inset: Midplane confocal image of fission yeast ex¬ 
pressing pomlp-tomato in the tip (red) and cdr2-GFP at 
the cell center (green). B Interactions between pomlp and 
cdr2p on the membrane. Pomlp actively disassociates cdr2p 
from the membrane (yellow box) and also inhibits the abil¬ 
ity of cdr2p to cluster (orange box). C Normalized concen¬ 
tration profiles for pomlp (red) and cdr 2 p (green) generated 
from Eqs. |H3| |12|13| Parameters for pomlp as Fig. For 
cdr 2 p, rji = 72/2 = Jc = 2 . 6 s“^, [P]o = 20pm~^, 

Dc,i = 0.5pm'^s~^, Dc ,2 = 0.01p,m?s~^, pc,i = Me, 2/2 = 
2.5 X 10“®s“^, Oc = O.lai and Pc = O.lji. D Positional error 
in the cdr2p profile when the pomlp insertion rate has 25% 
Gaussian variation. Results for clustering model (triangles) 
and equivalent SDD model (circles) shown. 


ior: ac(\P]) = Q:c,o(l + This simplified 

model recapitulates the observed pomlp and cdr2p pro¬ 
files, Fig. The fitting from the above model to the 
measured pomlp and cdr2p profiles represents a signifi¬ 
cant improvement over previous models [5S1 [321 [SD] both 
in replicating the spatial profiles and reproducing the dy¬ 
namics of clustering and gradient formation. 

To test the system robustness we created 200 pomlp 
profiles, with insertion rate normally distributed with 
standard deviation 25% of the mean and subsequent 
cdr2p profiles using Eqs. PPl At each position investi¬ 
gated we found the mean total cdr2p concentration and 
used this to define the threshold concentration for that 
position. For each individual cdr2p prohle we then mea¬ 
sured the position where it had each particular threshold 
concentration and hence calculated the positional preci¬ 
sion of cdr2p specification by calculating the standard de¬ 
viation in these positions. Near the source there is large 
variation due to the big intensity changes in pomlp be- 
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tween cells. Near the cell center there is large error due 
to the pomlp profile becoming increasingly flat. How¬ 
ever, around 2 — 3^m from the source which corresponds 
to the region where the boundary between pomlp and 
cdr2p is defined - we see that the cdr2p can be positioned 
more accurately (compared with equivalent SDD model) 
by pomlp when it clusters, Fig. Up. Of course, includ¬ 
ing intrinsic fluctuations would decrease the positional 
precision so the given accuracy represents a best case 
scenario (experimentally, errors of around l^m are typi¬ 
cally observed [15] )■ In conclusion, a dynamic clustering 
mechanism for concentration gradient formation can pro¬ 
vide robust positional information on relevant time and 
spatial scales for a biologically plausible scenario. 

V. DISCUSSION AND CONCLUSIONS 

Previous modelling of clustering within cells has pre¬ 
dominantly analyzed receptor clustering. Here, we have 
focused on the role of clustering in the formation of con¬ 
centration gradients and demonstrated that subcellular 
gradients can be formed via clustering on realistic spatial 
and temporal scales. This work represents a significant 
advance on previous models of subcellular concentration 
gradient formation [29l |46l EO] [52] as it accounts for pro¬ 
tein clustering and diffusion in a mechanistic (though still 
relatively straightforward) framework that also incorpo¬ 
rates realistic protein dynamics, and allows predictions to 
be made about the behaviour of specific dynamic com¬ 
ponents (see below). 

Our modeling enables the following predictions regard¬ 
ing the clustering of signaling molecules, (i) the process 
of clustering does not significantly stabilize (whether that 
be by extending protein lifetime or membrane association 
time) the individual molecules within the concentration 
gradient. If clustering significantly stabilizes the protein 
in the relevant system then the resultant chemical gradi¬ 
ent is not robust. This result is qualitatively consistent 
with the dynamics of pomlp, where the cluster lifetime 
is significantly shorter than the pomlp lifetime on the 
membrane |55|. (ii) The ratio of the fragmentation rate 
to disassociation rate plays an important role in the gra¬ 
dient formation. Either too small or big a ratio results in 
reduced spatial diffusion modulation and hence less ro¬ 
bust concentration gradients. Therefore, systems that 
use clustering in concentration gradient formation are 
likely to have carefully tuned fragmentation and disas¬ 
sociation rates and experimental perturbation of either 
should result in significantly reduced robustness, (hi) 
Clustering is favorable in systems that have a single de¬ 
cision to make (e.g. placement of division boundary) but 
it is less likely to be used in systems that specify multiple 
boundaries. Near the source there is increased inaccu¬ 
racy due to larger fluctuations in the concentration and 
at very large distances the profile becomes very flat due 
to only small, fast, clusters being present. Depending on 


the specific parameters, there will likely be an optimal 
region for concentration gradient interpretation. 

Given the advantages described above, why is cluster¬ 
ing of signaling molecules {i.e. morphogens) not observed 
more commonly in multicellular systems? In single cells 
where a simple decision is made by the signaling pathway 
{e.g. where to define the cell center) then clustering may 
be advantageous as dynamic parameters can be tuned 
to maximize precision at the relevant position but for 
morphogen gradients, that typically define three or more 
threshold positions across their spatial range, clustering 
may not be benehcial. Further, multicellular organisms 
typically have more time and complexity to adjust for 
variation in the input signal, such as via feedback net¬ 
works [T] ESI . We note that the Hedgehog signaling pro¬ 
tein is observed to cluster jSlj but this is likely due to its 
need for a chaperone to aid it in traversing through the 
intercellular space due to its hydrophobic nature [54] . 

We have considered a one-dimensional mean-field sce¬ 
nario. Two-dimensional simulations of the clustering 
model could be interesting, as the stochastic noise in such 
a (plausible) scenario is non-trivial, particularly near the 
insertion region. However, experimental evidence sug¬ 
gests that the approximations used here are relevant in 
the signaling region {i.e. away from the source). Fur¬ 
ther, extension to Smoluchowski processes where clusters 
directly interact with each other may be interesting but, 
as noted above, such a scenario is not consistent with 
current experimental observations of subcellular concen¬ 
tration gradients. 

Overall, we have presented a quantitative framework 
for understanding subcellular concentration gradient for¬ 
mation. In particular, our model can simultaneously 
replicate observed experimental spatial profiles and dy¬ 
namics. Importantly, despite the apparent complexity of 
such a clustering model, through biophysical arguments 
we reduced our model to two relevant parameters to de¬ 
scribe the concentration gradient. The resulting concen¬ 
tration profile is robust to relevant biochemical fluctua¬ 
tions so long as the process of clustering does not signif¬ 
icantly stabilize the signaling molecule. 
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